/*
 *  matrix.cpp
 *  T3nsors2
 *
 *  Created by Michael Barriault on 10-06-22.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "matrix.h"
#include "domain.h"
#include "scalar.h"
#include "vector.h"
#include "others.h"
#include "macros.h"
#include <iostream>
#include <cstdlib>

namespace t3 {
	
	matrix::matrix(void) {
		return;
	}
	
	matrix::~matrix(void) {
		FOR(n,N) T[n]->release();
		delete[] T;
		O->release();
		if ( id != "" )
			print((string)"Deleting " + id);
	}
	
	matrix* matrix::init(domain* O, int N, string id) {
		matrix* M;
		try {
			M = new matrix;
			//M->T = (t3::vector**)malloc(sizeof(vector*)*N);
			M->T = new vector*[N];
			FOR(n,N) M->T[n] = vector::init(O, N, id+itos(n));
		}
		catch (...) {
			Alloc_Error(id);
			throw;
		}
		M->O = O;
		M->O->retain();
		M->N = N;
		M->id = id;
		M->count = 1;
		print((string)"Making new" + id);
		return M;
	}
	
	matrix* matrix::copy(matrix* M) {
		if ( this->O == M->O ) FOR(n,this->N) {
			this->c(n)->copy(M->c(n));
		}
		return M;
	}
	
	matrix* matrix::retain(void) {
		print(id + "Going from " + itos(count) + " to ",true);
		count++;
		print(itos(count));
		return this;
	}
	
	void matrix::release(void) {
		print(id + "Going from " + itos(count) + " to ",true);
		count--;
		print(itos(count));
		if ( count <= 0 ) delete this;
	}
	
	vector*& matrix::c(int a) {
		return T[a];
	}
	
	scalar*& matrix::c(int a, int b) {
		return T[a]->c(b);
	}
	
	matrix* diff(vector* x) {
		matrix* dx = matrix::init(x->O, x->O->N, (string)"d"+x->id);
		FOR(p,dx->N) {
			dx->c(p)->release();
			dx->c(p) = diff(x->c(p));
		}
		return dx;
	}
	
	scalar* div(vector* x) {
		matrix* Dx = diff(x);
		scalar* dx = scalar::init(x->O, (string)"div"+x->id);
		FOR(a,x->N) FOR(o,x->O->size()) dx->at(o) += Dx->c(a,a)->at(o);
		Dx->release();
		return dx;
	}
	
	vector* curl(vector* x) {
		if (x->N != 3)
			return x;
		matrix* Dx = diff(x);
		vector* dx = vector::init(x->O, 3, (string)"curl"+x->id);
		FOR(o,x->O->size()) {
			dx->c(0)->at(o) = Dx->c(2,1)->at(o)-Dx->c(1,2)->at(o);
			dx->c(1)->at(o) = Dx->c(0,2)->at(o)-Dx->c(2,0)->at(o);
			dx->c(2)->at(o) = Dx->c(1,0)->at(o)-Dx->c(0,1)->at(o);
		}
		Dx->release();
		return dx;
	}
}
